/*******************************************************************************

	Blackman & Villalobos (2020) 
		USE FORESTS OR LOSE THEM? REGULATED TIMBER 
		EXTRACTION AND TREE COVER LOSS IN MEXICO  

	
	NOTES:
	1. This dofile replicates all tables and figures 
	2. Written by: Laura Villalobos
	3. Date: June 26th, 2020

-------------------------------------------------------------------------------*/


********************************************************************************
********************************************************************************

//set seed for replication purposes
	set seed 555
	set sortseed 555

//prepare data
	global root "C:\Users\ALLENB\Documents\A PROJECTS\Mexico permits and deforestation\PAPER\REPLICATION data and code"															// <------------------- Replace ".." with the repository's root folder
	use "$root\Blackman_Villalobos_JAERE_2020.dta", replace
	
* ------------------------------------------------------------------------------
* 								Define globals
* ------------------------------------------------------------------------------
	#delimit ;
	// all covariates;
		global covariates 																		
			"scaled_cropprice scaled_pdens scaled_indigenous scaled_ejido						
			phina_ejidatcom phina_avecindados phina_posesionarios phina_publicationyear			
			scaled_oppo scaled_marg_index00 													
			scaled_pa parcelized   																
			scaled_temp scaled_precp scaled_wc_tmean scaled_wc_prec 							
			scaled_dem scaled_aspect scaled_slope scaled_area_ha 								
			scaled_distocity scaled_distoclear scaled_dist_SEMARNAT 						
			scaled_carbon treecover 															
			share_area_ha_70s share_area_ha_1993 share_area_ha_2000 							
			biome_1 biome_2 biome_3 biome_13 biome_14";									

	// time invariant covariates;
		global timeinv 																				
			"scaled_ejido scaled_oppo scaled_marg_index00 											
			phina_ejidatcom phina_avecindados phina_posesionarios phina_publicationyear			
			scaled_pa parcelized  																
			scaled_wc_tmean scaled_wc_prec 														
			scaled_dem scaled_aspect scaled_slope scaled_area_ha 								
			scaled_distocity scaled_distoclear scaled_dist_SEMARNAT 						
			scaled_carbon treecover 															
			share_area_ha_70s share_area_ha_1993 share_area_ha_2000 							
			biome_1 biome_2 biome_3 biome_13 biome_14";						

	// time variant covariates;
		global timevar 																			
			"scaled_cropprice scaled_pdens scaled_indigenous scaled_temp scaled_precp ";	

	// socioeconomic covariates;
		global socioeconomic 																	
			"scaled_cropprice scaled_pdens scaled_indigenous scaled_ejido						
			phina_ejidatcom phina_avecindados phina_posesionarios phina_publicationyear				
			scaled_oppo scaled_marg_index00"; 
		
	// institutional covariates;
		global institutional 																	
			"scaled_pa parcelized"; 

	// climatological covariates;
		global climatological 																	
			"scaled_temp scaled_precp scaled_wc_tmean scaled_wc_prec"; 

		
	// geophysical covariates;
		global geophysical "scaled_dem scaled_aspect scaled_slope scaled_area_ha 				
			scaled_distocity scaled_distoclear scaled_dist_SEMARNAT 						
			scaled_carbon treecover 															
			share_area_ha_70s share_area_ha_1993 share_area_ha_2000 							
			biome_1 biome_2 biome_3 biome_13 biome_14";

	// region covariates;
		global region 																			
			"region_1 region_2 region_3 region_4 region_5";
			
	// time invariant covariates without pre-treatment deforestation vars;
		global timeinv_long_run 
			"scaled_ejido scaled_oppo scaled_marg_index00
			phina_ejidatcom phina_avecindados phina_posesionarios phina_publicationyear
			scaled_pa parcelized
			scaled_wc_tmean scaled_wc_prec
			scaled_dem scaled_aspect scaled_slope scaled_area_ha 
			scaled_distocity scaled_distoclear scaled_dist_SEMARNAT
			scaled_carbon treecover 
			biome_1 biome_2 biome_3 biome_13 biome_14";

					
	#delimit cr

* ------------------------------------------------------------------------------
* 								xtset 
* ------------------------------------------------------------------------------
	
	xtset fmu_id year	

* ------------------------------------------------------------------------------
* 					Propensity score matching by region
* ------------------------------------------------------------------------------
	
	gen u=uniform() 					   											// gen a random number

	foreach num of numlist 1 2 3 4 5  {
		sort u																		// sort randomly before using psmatch2
		psmatch2 evertreated $timeinv if tag==1 & region==`num', n(1) noreplacement	// match FMUs (cross-section)
		gen psm_`num'_area=1 if _weight!=. & _supp==1  								// identify the matched sample
		sort clv_unica year															// sort database by FMU year
		bysort clv_unica: carryforward psm_`num'_area, replace 						// copy weights (panel)
		}

	egen psm_append=																/// pool sample of matched observations per region
		rowmax(psm_1_area psm_2_area psm_3_area psm_4_area psm_5_area) 		
	

* ------------------------------------------------------------------------------
* 						TABLE 1: silvicultural system
* ------------------------------------------------------------------------------
	
	sum  sistemasilvicola_* if evertreated==1 & tag==1
		
* ------------------------------------------------------------------------------
* 						TABLE 3: summary statistics
* ------------------------------------------------------------------------------
	
	//part I. means and t-test

		foreach var in loss in_permit $socioeconomic $institutional $climatological $geophysical $region {
			sum `var'																		//col 1							
			sum `var' if evertreated==1														//col 2
			sum `var' if evertreated==0														//col 3
			sum `var' if psm_append==1 														//col 6
			
			sum `var' if psm_append==1 & evertreated==1										//col 7
			local mt=r(mean)																//to calculate bias before matching "by hand"
			local vt=r(Var)																	//to calculate bias before matching "by hand"
			
			sum `var' if psm_append==1 & evertreated==0										//col 8
			local mc=r(mean)																//to calculate bias before matching "by hand"
			local vc=r(Var)																	//to calculate bias before matching "by hand"
			
			ttest `var', unequal by(evertreated)											//col 5
			ttest `var' if psm_append==1 , unequal by(evertreated)							//col 10
			
			display ((`mc'-`mt') / sqrt((`vt'+`vc')/2)*100)	 								//col 4 bias before matching "by hand"
			pstest `var', t(evertreated) summary mweight(psm_append) support(psm_append)	//col 9
		}


	//part II. Median standardized bias  
	
		//unmatched
			pstest $covariates $region, t(evertreated) raw summary  
		
		//matched
			pstest $covariates $region, t(evertreated) summary  mweight(psm_append) support(psm_append) 

		
* ------------------------------------------------------------------------------
* 		TABLE 4: FMU-level cross-sectional propensity score probit results
* ------------------------------------------------------------------------------

	foreach num of numlist 1 2 3 4 5  {
		
		probit evertreated $timeinv if tag==1 & region==`num'

		margins if tag==1 & region==`num', dydx($timeinv) atmeans post

	}

	
* ------------------------------------------------------------------------------
* 						TABLE 5: main results
* ------------------------------------------------------------------------------

	//unmatched
		
		//FE, all regions 
			xtreg loss in_permit $timevar i.year, fe cluster(clv_unica)
			
		//MDE
			scalar MDE= _se[in_permit]*2.8
			display MDE
			
		//counterfactual
			margins, at(in_permit = 0) post
			
		//MDE/counterfactual
			display MDE/_b[_cons]

	//matched
		
		//FE, all regions matched
			xtreg loss in_permit $timevar i.year if psm_append==1, fe cluster(clv_unica)

		//MDE
			scalar MDE= _se[in_permit]*2.8
			display MDE
			
		//counterfactual
			margins, at(in_permit = 0) post
				
		//MDE/counterfactual
			display MDE/_b[_cons]
	
* ------------------------------------------------------------------------------
* 				TABLE 6: treatment effect heterogeneity by region 
* ------------------------------------------------------------------------------

	foreach num of numlist 1 2 3 4 5  {
		
		//main regression
			xtreg loss c.in_permit c.in_permit#i.region_`num' $timevar i.year if psm_append==1, fe cluster(clv_unica)
			est store Region_`num'
		
		//marginal effects 
			margins, dydx(in_permit) at(region_`num' = 1) post 						 
		
		//baseline probability 
			est restore Region_`num'
			margins, at(in_permit = 0) post											
	}

* ------------------------------------------------------------------------------
* 			TABLE 7: Treatment effect heterogeneity by FMU characteristic 
* ------------------------------------------------------------------------------

	foreach var of varlist unsscaled_oppo scaled_marg_index00 scaled_pa unsscaled_wc_prec biome_13 {
		
		sum `var'
		local mean=round(r(mean))												
		
		//main regression
			xtreg loss c.in_permit c.in_permit#c.`var' $timevar i.year if psm_append==1, fe cluster(clv_unica)
			est store cov_`var'
		
		//marginal effect evaluated at the mean
			margins, dydx(in_permit) at(`var'=`mean') post 							 
		
		//baseline probabiliy 
			est restore cov_`var'
			margins, at(in_permit = 0) post				   							
	}

		
********************************************************************************
********************************************************************************
								
						*SUPPLEMENTARY MATERIALS*
								
********************************************************************************
********************************************************************************
		
* ------------------------------------------------------------------------------
* 					TABLE A1: New permits by year and state  
* ------------------------------------------------------------------------------
	
	tab cve_ent year if new_permit==1, matcell(cell) matrow(rows)

* ------------------------------------------------------------------------------
* 			TABLE A2: overall, within-group, and between-group variation   
* ------------------------------------------------------------------------------

	xtsum loss in_permit
	xtsum loss in_permit if psm_append==1 

* ------------------------------------------------------------------------------
* 					TABLE A3: testing the common trends assumption 
* ------------------------------------------------------------------------------


	foreach num of numlist 1 2 3 4 5  {
		sort u
		psmatch2 evertreated $timeinv_long_run if tag==1 & region==`num', n(1) noreplacement
		gen psm_`num'_long_run=1 if _weight!=. & _supp==1
		tab _supp
		sort clv_unica year
		bysort clv_unica: carryforward psm_`num'_long_run, replace
	}

	egen psm_append_long_run=rowmax(psm_1_long_run psm_2_long_run psm_3_long_run psm_4_long_run psm_5_long_run) 
	sort clv_unica year
	bysort clv_unica: carryforward psm_append_long_run, replace
	
	//delta70_93
		reg delta70_93 evertreated if tag==1  											//unmatched
		reg delta70_93 evertreated $timeinv_long_run if tag==1 							//unmatched, controls
		reg delta70_93 evertreated if tag==1 & psm_append_long_run==1					//matched
		reg delta70_93 evertreated $timeinv_long_run if tag==1 & psm_append_long_run==1	//matched controls
	
	//delta93_00
		reg delta93_00 evertreated if tag==1
		reg delta93_00 evertreated $timeinv_long_run if tag==1 
		reg delta93_00 evertreated if tag==1 & psm_append_long_run==1
		reg delta93_00 evertreated $timeinv_long_run if tag==1 & psm_append_long_run==1
			
	//delta93_00
		reg delta70_00 evertreated if tag==1
		reg delta70_00 evertreated $timeinv_long_run if tag==1 
		reg delta70_00 evertreated if tag==1 & psm_append_long_run==1
		reg delta70_00 evertreated $timeinv_long_run if tag==1 & psm_append_long_run==1


* ------------------------------------------------------------------------------
* 					TABLE A4: correlation coeffients   
* ------------------------------------------------------------------------------
	
	correlate region_2 unsscaled_oppo scaled_marg_index00 unsscaled_wc_prec 

* ------------------------------------------------------------------------------
* 				TABLE A5: T-test of distribution of covariates 
* ------------------------------------------------------------------------------

	gen high_opp=1 if unsscaled_oppo>90
	replace high_opp=0 if unsscaled_oppo<=90

	gen high_marg=1 if scaled_marg_index00>3
	replace high_marg=0 if scaled_marg_index00<=3

	gen high_ppp=1 if unsscaled_wc_prec>110
	replace high_ppp=0 if unsscaled_wc_prec<=110


	foreach var in high_opp high_marg high_ppp{

		//means
			sum `var' if region_2==0
			sum `var' if region_2==1

		//t-test
			ttest `var', unequal by(region_2)

	}		
			
		
* ------------------------------------------------------------------------------
* 						TABLE A6: robustness checks
* ------------------------------------------------------------------------------
	
	//FIRST 5 YEARS
		xtreg loss first_5y $timevar i.year if psm_append==1, fe cluster(clv_unica)

	//FIRST 2 YEARS
		xtreg loss first_2y $timevar i.year if psm_append==1, fe cluster(clv_unica)

	//ANTICIPATION EFFECTS
		xtreg loss in_permit inpermit_t_2 $timevar i.year if psm_append==1, fe cluster(clv_unica)

	//CONTROL FOR SPILLOVER
		
		//PSM by region
			foreach num of numlist 1 2 3 4 5  {
				sort u
				psmatch2 evertreated $timeinv if tag==1 & region==`num' & drop_buffer5km!=1, n(1) noreplacement  
				gen psm_`num'_area_5km=1 if _weight!=. & _supp==1						
				sort clv_unica year														
				bysort clv_unica: carryforward psm_`num'_area_5km, replace
			}

		//pool sample of matched observations per region
			egen psm_append_5km=rowmax(psm_1_area_5km psm_2_area_5km psm_3_area_5km psm_4_area_5km psm_5_area_5km)  
			sort clv_unica year
			bysort clv_unica: carryforward psm_append_5km, replace

		//main regression
			xtreg loss in_permit $timevar i.year if psm_append_5km==1, fe cluster(clv_unica) 


	//TREATED ONLY SAMPLE 
		xtreg loss in_permit $timevar i.year if evertreated==1, fe cluster(clv_unica)

	//MATCHING WITH REPLACEMENT & N=3 
		foreach num of numlist 1 2 3 4 5  {
			sort u
			psmatch2 evertreated $timeinv if tag==1 & region==`num', n(3)
			gen psm_`num'_area_m3=1 if _weight!=. & _supp==1
			gen weight`num'=_weight
			sort clv_unica year
			by clv_unica: carryforward psm_`num'_area_m3 if region==`num', replace
			by clv_unica: carryforward weight`num' if region==`num', replace
		}

		//pool sample of matched observations per region
			egen psm_append_m3=rowmax(psm_1_area_m3 psm_2_area_m3 psm_3_area_m3 psm_4_area_m3 psm_5_area_m3) 
			egen weight_append_m3=rowmax(weight1 weight2 weight3 weight4 weight5) 
		
		//main regression
			xtreg loss in_permit $timevar i.year [pweight=weight_append_m3] if psm_append_m3==1, fe cluster(clv_unica)

			
	//DROP INSIGNIFICANT COVARIATES 	

		foreach num of numlist 1 2 3 4 5  {

			probit evertreated $timeinv if tag==1 & region==`num'

			margins if tag==1 & region==`num', dydx($timeinv) atmeans post
		}

		#delimit ;
			local timeinv_test_r1 " 
				scaled_wc_tmean 
				scaled_area_ha 
				scaled_dist_SEMARNAT 
				scaled_carbon 
				biome_1 biome_2 biome_3 biome_13 biome_14";

			local timeinv_test_r2 "scaled_ejido scaled_oppo  
				scaled_pa parcelized  
				scaled_wc_tmean 
				scaled_dem scaled_area_ha 
				scaled_distocity scaled_dist_SEMARNAT 
				scaled_carbon treecover 
				share_area_ha_70s 
				biome_1 biome_2 biome_3 biome_13 biome_14";

			local timeinv_test_r3 " 
				share_area_ha_70s";

			local timeinv_test_r4 "scaled_oppo  
				phina_ejidatcom phina_avecindados phina_posesionarios scaled_wc_prec
				scaled_distocity scaled_carbon";

			local timeinv_test_r5 "scaled_oppo  
				scaled_pa 
				scaled_wc_tmean 
				scaled_dem scaled_area_ha 
				scaled_distocity scaled_distoclear  
				scaled_carbon treecover 
				share_area_ha_2000 
				biome_1 biome_2 biome_3 biome_13 biome_14";
		#delimit cr

		foreach num of numlist 1 2 3 4 5  {
			sort u
			psmatch2 evertreated `timeinv_test_r`num'' if tag==1 & region==`num', n(1) noreplacement
			gen psm_`num'_area_test=1 if _weight!=. & _supp==1
			gsort clv_unica year
			bysort clv_unica: carryforward psm_`num'_area_test, replace
		}

		egen psm_append_test=rowmax(psm_1_area_test psm_2_area_test psm_3_area_test psm_4_area_test psm_5_area_test) 
		
		//main regression
			xtreg loss in_permit $timevar i.year if psm_append_test==1, fe cluster(clv_unica)
			
				
	//STEPWISE SPECIFICATION

		foreach num of numlist 1 2 3 4 5  {
			psestimate evertreated if tag==1 & region==`num', totry($timeinv) noquad  	 //select the variables, only linear
			sort u
			psmatch2 evertreated `r(K_l)' if tag==1 & region==`num', n(1) noreplacement  //Diferent set for each region
			gen psm_`num'_imbens=1 if _weight!=. & _supp==1								 //identify the matched sample
			sort clv_unica year
			by clv_unica: carryforward psm_`num'_imbens if region==`num', replace
		}

		//pool sample of matched observations per region
			egen psm_append_imbens=rowmax(psm_1_imbens psm_2_imbens psm_3_imbens psm_4_imbens psm_5_imbens) 
		
		//regression
			xtreg loss in_permit $timevar i.year if psm_append_imbens==1, fe cluster(clv_unica)
		
		//revise that No second order covariates are chosen!
			#delimit ;
			global imbensr1 "scaled_carbon biome_1 scaled_dist_SEMARNAT 
				scaled_area_ha scaled_wc_tmean phina_publicationyear scaled_oppo 
				share_area_ha_2000 scaled_distocity phina_avecindados parcelized";
			
			global imbensr2 "scaled_carbon biome_3 scaled_oppo share_area_ha_70s 
				scaled_pa scaled_wc_tmean treecover biome_1 scaled_area_ha 
				scaled_dist_SEMARNAT scaled_distocity parcelized scaled_ejido 
				scaled_aspect scaled_dem biome_13 share_area_ha_2000 
				scaled_marg_index00";
				
			global imbensr3 "share_area_ha_1993 scaled_oppo biome_2 parcelized 
				share_area_ha_70s scaled_area_ha biome_3";
				
			global imbensr4 "scaled_carbon phina_ejidatcom scaled_dist_SEMARNAT
				phina_avecindados scaled_wc_prec scaled_ejido biome_2 
				scaled_distocity share_area_ha_1993 share_area_ha_2000";
				
			global imbensr5 "treecover biome_3 scaled_area_ha share_area_ha_2000
				scaled_distoclear scaled_wc_tmean scaled_dem scaled_pa 
				scaled_distocity scaled_carbon phina_publicationyear 
				scaled_oppo scaled_dist_SEMARNAT scaled_aspect 
				phina_ejidatcom scaled_slope";
			#delimit cr

			psestimate evertreated if tag==1 & region==1, totry($imbensr1) nolin 
			psestimate evertreated if tag==1 & region==2, totry($imbensr2) nolin 
			psestimate evertreated if tag==1 & region==3, totry($imbensr3) nolin 
			psestimate evertreated if tag==1 & region==4, totry($imbensr4) nolin 
			psestimate evertreated if tag==1 & region==5, totry($imbensr5) nolin 

* ------------------------------------------------------------------------------
* 						TABLE A7: event study regression  
* ------------------------------------------------------------------------------
	
	xtreg loss permit_upto_3lags inpermit_lag_2 inpermit_lag0 - inpermit_lag5 permit_6pluslags ///
		$timevar i.year if psm_append==1, fe cluster(clv_unica) level(90) 
		est store event

	
	//Counterfactual
		margins, at(permit_upto_3lags=0 inpermit_lag_2=0 inpermit_lag0=0 ///
			inpermit_lag1=0 inpermit_lag2=0 inpermit_lag3=0 inpermit_lag4=0 inpermit_lag5=0 permit_6pluslags=0) post	
		
	//Three-year cumulative effects
		est restore event
		lincom _b[inpermit_lag0] + _b[inpermit_lag1] + _b[inpermit_lag2] 
		display r(estimate)
		display r(se)

		lincom _b[inpermit_lag3] + _b[inpermit_lag4] + _b[inpermit_lag5] 
		display r(estimate)
		display r(se)

	
********************************************************************************
********************************************************************************
								
								*FIGURES*
								
********************************************************************************
********************************************************************************



* ------------------------------------------------------------------------------
* 								FIGURE 4   
* ------------------------------------------------------------------------------

		preserve
			keep if tag==1 & psm_append_long_run==1
			collapse (mean) share_area_ha_70s share_area_ha_1993 share_area_ha_2000, by(evertreated)
			rename share_area_ha_70s share_area_ha_1970
			
			list
			reshape long share_area_ha_, i(evertreated) j(year)
			list
			
			graph set window fontface "Times New Roman"
			
			line share year if evertreated==1 ||   				 ///
			line share year if evertreated==0, lpattern(dash)    ///
			graphregion(fcolor(white) style(none)  				 ///
			color(gs16)) title("")  		     				 ///
			ytitle("Forest cover (%)") xtitle("")	     		 ///
			legend(order(1 "Ever-permitted" 2 "Never-permitted"))
		restore
		
* ------------------------------------------------------------------------------
* 								FIGURE 5   
* ------------------------------------------------------------------------------
	
	foreach var of varlist unsscaled_oppo unsscaled_wc_prec scaled_marg_index00  {

		local vtext : variable label `var'

		sum `var'
		local min=round(r(min))
		local sd=round(r(sd))
		local max=round(r(max))
		local mean=round(r(mean))

		//main regression
			xtreg loss c.in_permit c.in_permit#c.`var' $timevar i.year if psm_append==1, fe cluster(clv_unica) ///for margins to work, specify in_permit as continuous  

		//marginal effects
			margins, dydx(c.in_permit) at(`var'=(`min'(`sd')`max'))
		
		//marginal effects graphs
			marginsplot, 															///
				xlabel(`min'(`sd')`max', labsize(small)) xlabel(#6) 				///
				xtitle("`vtext'", size(small)) 										///
				ylabel(`mylab', labsize(small)) yline(0)							///
				ytitle("Marginal effect on forest loss", size(small)) 				///
				recast(line) recastci(rarea) 										///
				ci1opts(fintensity(30)) ciopt(color(%20))  							///
				graphregion(fcolor(white) style(none) color(gs16)) 					///
				title("") level(90)													

		//density graphs
			twoway 																		///
				(kdensity `var' if region_2==1,  										/// 
					color(red) 															///				
					xlabel(#6) xlabel(, labsize(small)) xtitle("`vtext'", size(small))	///
					ylabel(, labsize(small))  ytitle("Kernel density", size(small))) 	///
				(kdensity `var' if region_2==0,   										/// 
					color(blue) 														///
					xlabel(#6) xlabel(, labsize(small)) xtitle("`vtext'", size(small))	///
					ylabel(, labsize(small)) ytitle("Kernel density", size(small))) 	///
				(kdensity `var' ,    													///
					color(black) 														///
					xlabel(#6) xlabel(, labsize(small)) xtitle("`vtext'", size(small))	///
					ylabel(, labsize(small)) ytitle("Kernel density", size(small))) , 	///
				legend(order(1 "South region" 2 "Other regions" 3 "All FMUs") 			///
				size(small) width(100) col(3)) 											///
				graphregion(fcolor(white) style(none) color(gs16))
		}

			
* ------------------------------------------------------------------------------
* 								FIGURE A1   
* ------------------------------------------------------------------------------
			
	preserve
		
		est replay event, level(90)
		
		//obtain matrix with coefficients and CI
			matlist r(table)														
			mat coef_ci = r(table)
				
		//saves this matrix as variables in the dataset
			svmat  coef_ci	
				
		//keep only these numbers	
			keep coef_ci1 - coef_ci9												
			gen seq=_n
				
		//delete information in the matrix other than coeff and CI
			keep if seq==1 | seq==5 | seq==6		
				
		//gen an empty bin for the reference category
			gen coef_ci10 = .														
			
		//tidy up
			drop seq
			gen seq=_n
			gen stat="coef" if seq==1
			replace stat="cil" if seq==2
			replace stat="ciu" if seq==3
			
		//get it in the right shape for graphs
			reshape long coef_ci, i(stat) j(lag)									

		//re-organize information (note that the coef_ci are not labeled in the right sequence)
			replace lag=-3 if lag==1
			replace lag=-2 if lag==2
			replace lag=0 if lag==3
			replace lag=1 if lag==4
			replace lag=2 if lag==5
			replace lag=3 if lag==6
			replace lag=4 if lag==7
			replace lag=5 if lag==8
			replace lag=6 if lag==9
			replace lag=-1 if lag==10

		
			sort lag seq

		//graph
			graph set window fontface "Times New Roman"

			line coef_ci lag if seq==1, sort cmissing(n) lcolor(navy) 			|| 			///
			line coef_ci lag if seq==2, cmissing(n) lpattern(dash) lcolor(navy) || 			///
			line coef_ci lag if seq==3, cmissing(n) lpattern(dash) lcolor(navy) 	 		///
			xlabel(-3 "<=-3" -2 "-2"  -1 "-1"  0 "0" 1 "1" 2 "2" 3 "3" 4 "4" 5 "5" 6 ">=6") ///
			yline(0) ytitle("Regression coefficient")   									///
			xline(0) xtitle("Years since the permit was granted") 							///
			legend(order(1 "Estimated treatment effect" 2 "90% confidence interval")) 		///
			graphregion(fcolor(white) style(none) color(gs16)) title("")

	restore		
	
 
 log close	
	*****END
